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Abstract 

The  transport  of  particles  in  a  Particle-In-Cell  (PIC)  method  is  traditionally  handled  by  a 
staggered  algorithm,  second-order  accurate  in  time,  originally  developed  by  Boris  [1-2],  The 
scheme  is  very  efficient  and  although  it  is  stable  for  time  steps  large  compared  to  the  cyclotron 
period  (“gyro-period”),  it  ceases  to  be  accurate  in  that  case.  In  cases  of  strong  applied  magnetic 
field,  this  can  impose  an  impractical  time-step  restriction.  An  alternative  approach  is  to  average 
over  the  orbital  motion  and  consider  only  that  of  the  guiding-center;  this  has  led  to  so-called 
gyrokinetic  simulations  [3].  However,  that  approach  can  also  lead  to  some  inaccuracies,  due  to 
the  loss  of  information  regarding  the  phase  of  the  orbital  motion.  Furthermore,  it  may  also  be 
desirable  to  have  an  algorithm  that  is  not  staggered  in  time,  in  order  to  guarantee  exact 
conservation  of  total  energy  at  all  times.  In  this  paper,  we  present  an  algorithm  that  solves  the 
non-relativistic  equation  of  motion  exactly,  and  can  yield  exact  conservation  of  energy  for  large 
time  steps  (compared  to  gyroperid).  The  algorithm  accuracy  is  demonstrated  and  compared  with 
the  Boris  scheme.  These  preliminary  results  are  valid  for  the  homogenous  case  only,  and 
extension  to  spatially-varying  fields  should  be  considered  next. 
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1.  Introduction 

We  consider  the  problem  of  solving  the  non-relativistic  dynamical  equation  for  charged  particles 
in  arbitrary  electric  and  magnetic  fields: 

mv  -q  {E  +  vxB)  (1) 

This  is  the  basic  transport  process  in  Particle-In-Cell  (PIC)  codes,  which  is  usually  solved  using 
the  Boris  algorithm  [1],  defined  in  Appendix  A.  The  Boris  algorithm  is  a  computationally 
efficient  (i.e.  uses  a  minimum  number  of  operations)  algorithm,  second-order  accurate  in  time. 
Since  it  is  a  leap-frog  integrator,  it  is  also  usually  described  as  a  symplectic  algorithm,  i.e.  which 
conserves  a  discrete  analog  of  the  Hamiltonian  up  to  second-order  accuracy.  This  is  a  critically 
important  property  for  PIC  simulations,  which  usually  do  not  have  conservation  properties 
embedded  in  the  mathematical  formulation  as  in  continuum  models,  such  as  finite-volume  or 
finite-difference  schemes.  However,  it  is  important  to  exercise  some  caution  when  speaking  of 
energy  conservation  in  the  Boris  scheme;  as  a  leap-frog  algorithm,  it  uses  position  and  velocity 
staggered  in  time,  the  kinetic  and  potential  energies  are  not  computed  at  the  same  time.  After 
advancing  the  particle,  the  kinetic  energy  can  be  evaluated  from  the  velocity  field  at  time 
(n+P2) ,  while  the  potential  energy  can  be  obtained  exactly  from  the  particle  position  at  time 
(n+1) .  Thus,  the  kinetic  and  potential  energies  are  not  strictly  conserved  at  the  same  time. 
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Figure  1:  Schematic  of  leap-frog  Boris  algorithm. 


In  the  leap-frog  algorithm,  the  fields  are  used  at  the  mid-point  for  advancing  the  velocity,  i.e. 
fields  evaluated  at  time  (n)  are  needed  for  updating  the  velocity  from  to  .  In  an 

electrostatic  simulation,  the  electric  field  can  be  obtained  from  solving  Poisson’s  equation: 

VV  =  — (2a) 
^0 

and  (2b) 

In  (2a),  the  particle  density  at  a  given  location  (grid-point)  is  obtained  as  the  statistical  average  of 
the  contribution  of  neighboring  particle;  this  “scatter”  operation  maps  the  particles  onto  the  grid, 
and  various  interpolation  schemes  can  be  used  for  this  operation.  We  point  out  that  this  mapping 
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uses  the  partiele  loeations  at  time  level  (n) ,  and  therefore  the  eleetro-statie  potential  (and  eleetrie 
field)  are  naturally  synehronized  with  the  partiele  positions.  In  eleetro-statie  simulations  the 
magnetie  field  is  eonstant  and  there  is  no  eoneem  over  its  synehronization.  In  eleetro-magnetie 
simulations,  however,  both  fields  are  advaneed  in  time  and  the  proeedure  must  be  eonsistent  with 
the  Maxwell  equations: 


f  oV  •  F  =  -Pp 

(3a) 

BE  WxB  - 

"“a"-;.. 

(3b) 

—  =  -VxF 
dt 

(3c) 

V-B^O 

(3d) 

Equation  (3a),  where  is  the  eharge  density  from  the  partieles,  is  simply  Poisson’s  equation 
(2a);  in  (3b),  j ^  is  the  eurrent  density  from  the  partieles,  and  is  obtained  by  a  similar  mapping  of 
the  partiele  veloeities  onto  grid  points.  The  Maxwell  equations  are  naturally  synehronized  to 
seeond-order  aeeuraey  for  P*"\E*"*  and  jjQ^gygj-  ^j^g  leap-frog  algorithm  for 

particle  transport  is  of  the  form  =  f  ,  and  one  needs  to  interpolate  in 

time  one  of  the  fields  for  the  particle  push,  i.e.  B .  Note  that  the  leap-frog  algorithm  of  Figure  1  is 
not  the  unique  solution:  one  could  just  as  well  decide  to  choose  the  fields  r^"\B^"^  and 
-(«+V2)  ^^(«+v2)  Q^j^gj.  combinations  and  rely  on  the  time-interpolation  of  another  field  (matter 
or  particle)  to  re-establish  second-order  time  accuracy.  Higher-order  schemes  can  of  course  be 
obtained  with  iterative  methods. 


Figure  2:  Potential  positional  error  of  drift  dynamies  versus  gyro-radius 


The  Boris  algorithm  is  stable  at  high  values  of  the  magnetic  field,  i.e.  when  »  1 ,  although 
accuracy  is  lost  for  large  time  steps.  Practically  speaking,  the  time  step  in  PIC  simulations  using 
this  algorithm  is  restricted  such  that  tw^At « 1 ;  this  makes  the  scheme  highly  inefficient  in  cases 
of  strongly  magnetized  plasmas.  One  could  consider  an  alternative  approach  in  that  case,  where 
only  the  motion  of  the  guiding  center  is  modeled;  the  rotation  around  the  field  line  is  not  tracked, 
but  averaged  over  several  orbits.  This  “drift  dynamics”  approach  is  valid  when  the  gyro-radius 
the  characteristic  cell-size;  however,  significant  errors  can  be  introduced  even  when 
<  A^ .  Since  the  phase  of  the  gyro-motion  is  not  known  in  this  approximation,  the  particle 
position  is  effectively  randomized  on  a  scale  comparable  to  the  cell  size  (see  Figure  2).  This  can 
lead  to  errors  at  the  crossing  into  different  cells  or  boundaries,  and  errors  when  the  field  gradients 
on  the  scale  of  a  cell  size  are  non-negligible.  Therefore,  it  is  worth  investigating  the  construction 
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of  an  accurate  particle -push  algorithm  that  is  more  efficient  at  high  cyclotron  frequency,  yet 
remains  accurate  and  conserves  energy  to  a  high  level  of  accuracy.  This  is  the  object  of  the 
following  study. 

2.  Exact  Solution 

We  consider  here  the  case  of  constant  and  uniform  fields.  This  considerably  simplifies  the 
analysis  and  allows  us  to  obtain  an  exact  analytical  solution  to  the  non-relativistic  equations  of 
motion.  The  constant  field  approximation  is  valid  when  the  time-variation  of  the  fields  is 
neglected  during  the  time-step  (i.e.  first-order  time-accuracy  of  the  field  evolution);  the  extension 
to  higher-order  time-dependency  and  non-uniform  fields  will  be  examined  in  the  future.  We  will 
also  be  performing  a  transformation  to  the  reference  frame  aligned  with  the  magnetic  field.  Let  us 
first  define  the  laboratory  frame  (L)  by  the  italicized  letters  (x,y,z)  and  a  rotated  coordinate 

frame  by  fj,  such  that  the  unit  vector  is  aligned  with  the  magnetic  field,  i.e.  ^  =  b . 


Figure  3:  Reference  frame  transformation:  aligned  with  B  . 


The  rotation  operators  between  the  two  reference  frames  are  given  by: 


R 


^.V 

^z^ 

fc 

ScpCg 

\ 

-^e 

lx 

Iz 

0 

(4a) 

^z  ) 

A 

S<pSg 

rix 

4" 

A 

K 

c^s 

\ 

9 

= 

4'. 

■' 

S,pCg 

(4b) 

4'J 

0 

Cg 

) 

R  '  = 


where  we  have  used  the  condensed  notation  of  =  cos((9) ,  Sg  -  sin(0) ,  etc.  Once  the  rotation 
into  the  aligned  frame  is  performed  and  no  confusion  is  possible,  we  can  use  the  script  letters 
(x,  y,z)  to  denote  the  components  in  that  frame,  i.e.  (x,  y,  z)  =  (i^, .  We  will  also  denote 
vectors  in  that  frame  by  bold-face  type,  i.e.:  E  =  R  (H)  .E  . 


In  this  rotated  frame,  the  equation  of  motion  (1)  can  be  expressed  by: 
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(5) 


Vx  + — V 

m  m 

q  ^  qB 

Vy  =— E - 

m  m 

V,  =  — E, 


In  this  frame,  the  magnetic  field  has  a  component  only  in  the  z  direction,  and  therefore  =  5  , 
the  magnitude  of  the  magnetic  field.  We  can  define  a  normalized  electric  acceleration  field 
a  =  qEI m  and  the  cyclotron  frequency  co-qB /m  ,  which  is  a  signed  quantity.  The  solution  for 
the  z  -component  of  the  velocity  is  trivial  and  can  be  ignored  for  the  moment.  The  system  (5)  can 
be  reduced,  for  the  transverse  components,  to: 

^  =  ^  +  o]'5^  «’) 

where  the  underline  indicates  a  vector  in  the  transverse  directions  only.  An  additional  time 
derivative  of  (6)  yields  the  following: 

_  — 2,,  ,  (  0 


v  =  Y+  0  -a 


The  general  solution  of  (6-7)  is: 

Jvj,  sm{(M+(p)  +  {2L^I(o)sm{(M+(p)  +  {a ^1  Q)){\-cos{(M+(p)) 

-  1 V 0 cos{cot+(p)  -  (a  ! (o){\-cos{Q)t+(p))  +  (a  ^Ico)  sm{cot+(p)  ^ 

Let  us  denote  +  ^ .  It  is  to  verify  that  the  solution  (8)  satisfies  the  equations  of  motion: 


cos(<^)  +  a^  cos(<^)  +  aySin(<^)  ^  ((UVy +a^  q  co 

-  1-wv„sin(^^)-a^sin(^^)  +  a  cos(^^)  \-(aw ^+a  -  0 


Vq  sin(^^) -a^(ysin(^^)+a  (ycos(^^)  -co^v^+coa 


-CO  \ ^cos{(f>)  -a^(oco?,{(p)-a  (osm{(f>)  -co  \  -coa 


-  I  -  «  0  i  - 


0  (9a) 


0  (9b) 


Let  us  now  compute  the  solution  at  an  advanced  time  t  +  dt .  Lrom  (8)  we  have: 

fvg  sin(^^(i^^)  +  (ajco)  sin(^^(i^^)  +  {a ylco){y-cos{(jH-d(l))) 
Y(A  0  |v^cos(^^(i<^)  -  {ajco){\-cos{i^d4>))  +  {^yl(o)  sm{^d4>) 

Expanding  the  trigonometric  functions  we  find: 


y(f+£/f)  = 


Cg- K-ay/®  +  Sg- |vy+a^/w  +  {a^lQ})  ^ 

+cg-hy+^Jo^  -(ax/^y) 


O(^0-v(0 +  -[/", 

—  CO  Cg  V  Sg 


where  S  =  codt  and  ^{dt)  is  the  counter-rotation  matrix  around  the  magnetic  field: 

The  last  matrix  in  (1 1)  can  also  be  written  in  terms  of  this  rotation  matrix.  Let  us  define: 

fo  n  .  -1  -ifo  -\] 

o  =  (y|_j  qI  and  a  -co  A  ^ 

then  (11)  becomes: 


y(f+At)  =  fl(dt)  ■  y(f)  -t  ^ 

~^g  ^g~^ 


In  a  compact  form: 
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Av  =  v(?+A?)-v(?)  =  “  l]’ [y(0+'^  '-a]  (15) 

One  should  now  consider  the  case  of  vanishing  magnetic  field.  The  matrix  on  the  RHS  of 
(15)  is  singular  when  B  -  0 .  However,  it  can  be  combined  with  the  term  in  brackets  as  follows: 


[a(£/0-i]-o^'  = 


(  sin  A 

1-cos  A  ^ 

1  _ 

S 

S 

cos  5  -1 

sin  A 

1  A 

'■ - V" 

A  J 

•  dt 


(16) 


The  matrix  Aj  is  regular,  since: 

4i  ^ 


1  0 
0  1 


when  >0 


One  can  expand  (15)  to  the  next  order  in  A  =  ®At ,  leading  to: 


Av  = 


a ^dt+Q) dt{w ^+j?L ^dt)  -  \  {q) dt)^  \ ^ 


a  dt -CO dti\^+\di^dt)  -\{(odty\ 


(17) 


(18) 


We  see  that  the  second-order  accurate  {o{dt^))  solution  is  obtained  by  rotating  a  half-step 
advanced  solution,  as  expected.  We  can  also  verify  that  this  solution  is  identical  to  the  Boris 
algorithm,  by  comparing  (18)  with  (A.23b)  of  Appendix  A. 

The  opposite  limit  of  large  time  steps  compared  to  the  gyro-motion,  i.e.  codt  ^<x>,  is  also  of 
principal  interest.  In  that  case,  the  trigonometric  functions  oscillate  rapidly,  but  the  trajectory 
remains  bound.  One  can  perform  an  averaging  over  a  large  number  of  gyro-motions,  and 
eliminate  all  terms  proportional  to  these  functions  (<  cos  >=<  sin  >=  0 ).  The  remainder  is: 

0 


<  v{t+At)  >- 


(19) 


-Mco  0 

which  is  independent  of  the  time  step  At .  This  is  a  constant  velocity,  which  can  be  easily 
recognized  as  the  ExB  drift  velocity,  since  (19)  is  equivalent  to: 

ExB 

[-EJB 

Therefore,  the  formulation  (15)  automatically  recovers  the  drift  motion  of  the  guiding  center 
when  the  gyro-motion  is  not  resolved  -  with  a  randomized  rotation  around  the  magnetic  field. 

Let  us  now  look  at  the  exact  solution  for  the  particle  position.  From  (8),  we  obtain: 


<  v>= 


(20) 


1  1  1 

^-cos<^ 

-sin^^^ 

1  1 
•  3.  H - 

r  vi 

Q)  i 

.  sin<^  ^ 

0) 

^  sin^^ 

co' 

\~^xK 

+  Ao 


(21) 


The  expression  at  a  later  time  t  +  dt  can  be  expressed  as  function  of  the  original  phase  ^^(f)  and 
the  phase  difference  d(j)  =  5  -codt  by: 


x{t+dt)  =  — 

CO 


^  a^{t+dt) 
^-a^{t+dt)^ 


+  Ao  + 


-sj 


-(Vo/®)c^  -(a,/®^)c^  -{ay/(o^)s^ 
+(Vo/®)5^  -{a^lco^y^ 


(22) 


One  recognizes  again  the  rotation  matrix  (12)  in  that  expression,  which  allows  us  to  write  the 
displacement  as: 

^  a Jt\  ^^\-{-wJco)c^-{aJco^)c^-{aJco^)s^ 

+(Vo/®)s^  -{ajco^y^ 


Ax  =  l 

CO 


-ayt 


\^{dt)-\- 


(23) 


However,  from  (15)  one  can  also  recognize  the  following  expression: 
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\  ^-a^Jco 

y 

v,.+a  Z® 

>  ^ 


=  («-!)• 


+{aJco)s^  -{aylco)c^ 
VoC^  +(a,/®)c^  +{^ylo))s^ 


Inserting  (24)  into  (23),  we  finally  obtain: 


Ax  =  i. 

CO 


a^At -Av^ 


-a^At  +Av, 

For  the  displacement  along  the  magnetic  field,  the  exact  solution  is  of  course: 

AX||  =  V||(t)  At  +  yanAt^  =  V||(t+Y  At)  At 


(24) 

(25) 

(26) 


It  would  appear  that  the  transverse  displacement  (25)  has  a  singular  behavior  at  vanishing 
magnetic  field  strength,  due  to  the  «  '  factor.  However,  a  simple  Taylor  expansion  can  confirm 
that  this  is  not  the  case:  when  ®  — >•  0  one  can  use  (18)  into  (25)  to  verify  that,  as  expected: 


Ax; 


v^At  +  ja^At" 
v,.At  +  ^a„At^ 


+  o((y  ) 


(27) 


The  expression  (25)  is  therefore  valid  for  all  non-zero  values  of  the  magnetic  field.  We  can,  as 
before,  regularize  this  expression  in  the  case  of  a>-0;  after  some  simple  algebra,  we  obtain: 


Ax  = 


fl-c. 

5 

C<5-1 

5 

Ss 

•vAt  + 

5^ 

Sg-S 

5^ 

1  A 

S  ) 

1  A' 

) 

•  a  At^ 


(28) 


Ai  Aj 

We  have  recognized  the  first  matrix  (16),  and  defined  a  second  regularized  matrix  A2 .  Both  are 
finite  when  w  — >■  0 ,  since  in  that  case  (defining  also  the  following  S,  C  coefficients): 

^  S 
~6 


^0  - 


S 

‘  S  6 


^2-- 


l-Cx 


C2^- 
2  ^ 


1-Cx 


2 


S  2  "5^ 

Therefore,  the  procedure  outlined  above  is  applicable  in  all  cases  of  magnetic  field  values. 


(29a) 

(29b) 


3.  General  Algorithm 

One  can  construct  two  types  of  algorithms.  The  first  case  is  valid  only  for  5  >  s ,  i.e.  does  not 
require  regularization,  and  is  governed  by  the  following  operations: 

(1)  Transform  the  velocity,  position  and  acceleration  vectors  from  the  original  reference 
frame  into  the  rotated  frame  with  the  z  axis  aligned  with  the  magnetic  field. 

(2)  Compute  the  changes  to  the  transformed  velocity  vector,  separating  the  transverse  and 
parallel  components: 


Ay  =  (a((Z0-l)- 


v^-aylco 

v,,+a^/(y 


(30a) 
(30b) 

(3)  Compute  the  changes  to  the  transformed  position,  separating  the  transverse  and  parallel 
components: 


AV||  =  a||At 


Ax  =  i. 

CO 


ayM 

-a  At  +Av, 


(31a) 
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(31b) 


AX||  =  V||(?)  A?  +  ya||A?^ 

(4)  Transform  the  changes  back  into  the  original  frame  and  add  to  the  initial  values. 

In  the  second  case,  the  regularized  matrices  are  used  so  that  the  algorithm  remains  valid  for  all 
cases  of  field  values,  including  5  <s  .  Combining  transverse  and  parallel  components  ,  we  can 
express  the  velocity  change  as: 


(32a) 


(  -C 

0^ 

Cl 

0) 

a^At 

Av  = 

-^0 

-Co  0 

Vv 

+ 

-C, 

51 

0 

a  At 

1 0 

0 

oj 

Az. 

V  0 

0 

ij 

1 - 

N 

> 

A, 

C, 

0) 

\,a{ 

(c. 

-^2 

0 

1 

> 

Ax  = 

0 

WyAt 

+ 

S2 

C2 

0 

a^At 

V  0 

0 

ij 

v,At 

lo 

0 

a, At 

(32b) 


The  transformation  steps  (1)  and  (4)  remain  the  same. 

It  would  appear  that  the  proposed  scheme  is  very  expensive,  since  it  requires  the  evaluation  of 
several  matrices.  However,  in  the  case  of  constant  fields  studied  so  far,  these  matrices  can  be 
determined  once  the  fields  and  time  step  are  known.  The  transformation  matrices  (4)  can  be 
incorporated  into  the  definition  of  the  regularized  push  matrices  Aq,  Aj,  A2 ,  leading  to: 

Dk  =  -Ak  -R,  k  =  1,2,3  (33) 

The  changes  can  then  be  computed  directly  in  the  initial  (non-rotated)  frame: 


^  flAt 

Az;  =  Dj-z7  +  -2 — Dj-ii 

m 


and 


Ar  -  D^-vAt  + 


c\Ar 


m 


Dj  £ 


(34a) 

(34b) 


4.  Computational  Tests 

The  first  test  conducted  concerns  the  movement  of  a  single  particle  (a  positron)  in  static  fields; 
the  initial  velocity  is  null,  the  electric  field  is  1  kV/m  in  the  positive  y  -direction  and  the  magnetic 
field  is  1  Tesla  in  the  positive  z  -direction.  Under  such  initial  conditions  the  particle  executes  a 
cycloidic  movement  of  height  equal  to  where  is  the  drift  velocity 

in  the  x -direction.  The  motion  is  computed  for  three  cases  of  constant  time  steps,  being 
respectively  0.1/w^,  M  co^  and  10/®^.  Since  the  Boris  algorithm  requires  the  velocity  at  a  prior 

half-time  step,  that  initial  value  ( )  is  computed  from  the  exact  solution.  The  trajectories  for 
the  exact  solution,  the  Boris  algorithm,  and  the  regularized  algorithm  of  eqs.  (34)  are  shown  in 
Figure  4.  All  methods  are  in  very  good  agreement  for  small  time  step.  For  At-Mco^,  the  Boris 
algorithm  starts  to  show  some  noticeable  deviations  from  the  exact  solution;  first,  the  Larmor 
radius,  or  height  of  the  cycloid,  is  noticeably  larger;  second,  the  effective  gyro-frequency  is 
somewhat  lower,  leading  to  a  growing  de-phasing  with  the  exact  solution.  At  larger  time  steps 
(Figure  4c),  the  solution  from  the  Boris  algorithm  is  in  error  by  close  to  an  order  of  magnitude^. 
By  contrast,  the  current  algorithm  provides  a  solution  that  is  in  perfect  agreement  with  the  exact 
solution,  both  in  amplitude  and  phase. 


The  underline  is  eliminated,  sinee  these  are  now  3 -dimensional  variables. 

^  Note  that  the  magnitude  of  this  error  is  bounded,  i.e.  does  not  grow  in  time,  a  result  of  sympleetieity. 
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Y  [m]  Y  [m]  Y  [m] 


Figure  4:  Particle  trajectory  for  three  time  steps;  exact  solution  is  compared  to  the  results 
from  the  Boris  algorithm  and  the  current  scheme. 
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It  should  also  be  pointed  out  that,  despite  the  obvious  error  in  transverse  position  and  rotation 
frequeney  present  in  the  Boris  algorithm,  the  drifting  motion  in  the  x-direetion  is  accurately 
maintained.  This  is  evident  in  Figure  5,  where  the  difference  between  the  actual  position  and  the 
expected  position  from  the  constant  drift  is  plotted  versus  time.  The  natural  oscillations  are  due  to 
the  cycloid  motion  itself  and  of  amplitude  equal  to  the  Larmor  radius;  only  at  the  largest  time  step 
does  the  Boris  algorithm  deviates  from  the  expected  behavior.  The  current  method  (Figure  5b) 
yields  the  correct  drift  dynamics  at  all  time  steps. 


60x10'® - 1 - 1 - 1 - 1 - 

This  method 

>“  40-  - At  =  0.1/co^ 

^  -  At  =  1  /  (D^ 

-o  -  At=  10/(0^ 

20  - 

0  

-20-^ - ^ ^ ^ ^ - \- 

0.0  0.5  1.0  1.5  ,  ,  2.0  2.5x10'® 

time  [s] 

Figure  5:  X-position  versus  time,  normalized  to  theoretical  position  of  guiding  center  (Vnt). 

One  can  now  examine  the  impact  of  the  errors  on  conservation  properties,  i.e.  kinetic,  potential 
and  total  energies.  In  this  simple  test  case  with  an  imposed  external  field,  the  potential  energy  per 
mass  is  simply  -  {qElm)  ■  y ,  where  y  is  the  particle  position  along  the  y  -axis.  It  is 

important  to  point  out  that  for  the  Boris  algorithm,  the  kinetic  and  potential  energies  are  evaluated 
at  different  times,  i.e.: 

e%]  =  iqE/m)  y<«)  and  f  / 2  (35) 

Therefore  to  evaluate  the  total  energy  at  a  specific  time  (e.g.  one  must  interpolate  one  of  the 
variables  to  that  time  level.  Both  kinetic  and  potential  energies  are  shown  as  function  of  time  for 
the  Boris  algorithm  in  Figure  6. 


This  method 

- At  =  0.1/co^ 

-  At  =  1  /  (D^ 

-o  -  At=  10/(0^ 
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-3x10 


"T" 

1.2 


- r 

1.3 

X  [m] 

_ 1_ 


1.0 


1.1 


1.4x10'' 


^kin  Epot  Boris  (At  =  1/(dJ 


-3x10 


1.0 


"T" 

1.1 


1.2 


1.3 


1.4x10 


■9 


X  [m] 


Figure  6:  Kinetic  and  potential  energy  versus  distance  for  Boris  algorithm  -  all  cases  of  time 
steps.  Dashed  horizontal  lines  indicate  theoretical  limits  of  variation. 


It  can  be  seen  that  there  is  a  rapid  degradation  of  the  energy  conservation  as  the  time  step  is 
increased  to  values  of  the  same  order  or  beyond  the  gyro-period;  the  error  in  amplitude  of  the 
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particle  trajectory  leads  to  errors  in  potential  energy  which  become  severe  for  large  time  steps. 
Since  the  current  method  is  in  perfect  agreement  with  the  exact  solution  at  all  time  steps,  the 
energy  is  perfectly  conserved  in  that  case  (see  Figure  7). 


Figure  7:  Kinetic  and  potential  energy  versus  distance  for  current  method  -  all  time  steps  shown. 


To  evaluate  the  total  energy,  one  must 
account  for  the  dephasing  of  the 
velocity  and  position  in  the  case  of  the 
Boris  algorithm,  as  mentioned 
previously.  This  dephasing  can  be 
clearly  seen  when  plotting  both 
energies  versus  a  single  time  coordinate 
in  Figure  8.  The  shift  of  the  two 
curves  is  a  result  of  the  leap-frog 
algorithm.  One  can  correct  for  this  by 
plotting  the  kinetic  energy  versus  the 
proper  time  of  evaluation,  i.e.  the  set 
|^(«+V2)|  ^  Figure  7;  this  shifts 

all  the  points  to  the  left,  as  indicated  by 
the  black  arrows  of  Figure  8.  The  total 
energy  can  be  evaluated  at  the  set  of 
times  by  adding  the  potential 

energy  at  that  time  with  the  average  of 


Figure  8:  Kinetic  (red)  and  potential  (blue)  energies 
versus  a  single  time  coordinate. 


the  kinetic  energies  at  times  (?*"  v2)  ^^(n+i/2)^  interpolation  is  accurate  as  long  as  the  time 
step  is  sufficiently  small,  i.e.  such  a  linear  approximation  of  the  trajectory  between  the  two  times 
is  reasonable;  however,  it  can  lead  to  severe  errors  for  large  time  steps.  This  phase  error  of  the 
Boris  scheme  is  in  addition  to  the  amplitude  error  (effective  Larmor  radius)  already  observed,  for 
example,  in  Figure  6.  The  average  errors  on  total  energy  and  position  can  be  obtained  for  various 
values  of  the  time  step  and  magnetic  field.  The  error  on  the  total  energy  can  be  defined  here  as 
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err{E^^^)  =  {Ei^i-Ei^)/EQ,  where  is  the  initial  total  energy  (at  t  =  0)  and  Eq  is  a 
representative  energy  scale;  here,  Eq  -v\I2,  where  v^,  is  the  drift  velocity.  Similarly  for  the 
position,  the  error  is  defined  as  err(X)  =  s|x^-x^™‘"^|/Xo;  all  components  of  the  position  are 

contributing,  and  Xq  is  a  representative  length  scale  -  in  the  case  of  cold  particles  here,  the 
Larmor  radius,  i.e.  half  the  height  of  the  cycloid  motion.  Both  errors  are  shown  in  Figure  9  for  the 
two  schemes. 


Figure  9:  Error  on  energy  (a)  and  position  (b)  as  funetion  of  time  step  for  several  values  of  the 
magnetie  field,  for  the  Boris  and  our  algorithms. 
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Note  the  change  of  scale  between  the  left  (Boris  algorithm)  and  right  axis  (current  algorithm); 
clearly,  the  present  scheme  is  more  accurate  by  several  orders  of  magnitude. 

5.  Conclusions 

We  have  successfully  implemented  and  tested  a  new  particle  pusher  algorithm  that  can  effectively 
be  used  for  large  time  steps,  much  larger  than  the  gyroperiod;  the  method  is  based  on  the  exact 
solution  of  the  equations  of  motion,  but  does  not  require  tracking  the  phase  of  the  particle  motion 
around  the  field  line.  The  method  can  be  applied  to  arbitrarily  large  time  steps  and  yields  exact 
(down  to  machine  accuracy)  conservation  of  energy  and  exact  position.  The  method  is  currently 
restricted  to  the  non-relativistic  case,  and  to  uniform  fields.  Extension  to  the  relativistic  regime 
would  be  very  difficult,  since  there  is  no  longer  an  analytical  solution;  extension  to  the  non- 
uniform  (magnetic)  field  does  not  present  a-priori  any  difficulties,  but  this  must  be  verified. 

It  would  a-priori  appear  that  the  algorithm  is  computationally  expensive,  but  this  is  not 
necessarily  the  case.  The  push  matrices  (eqs.  32)  need  to  be  computed  only  once  for  each  time- 
step,  but  are  the  same  for  each  particle  in  this  case  of  uniform  field.  Thus,  the  method  would  be 
efficient  when  computing  a  large  number  of  particles  in  such  configurations,  e.g.  Penning  traps. 
In  the  case  of  weakly  non-uniform  fields,  one  can  also  attempt  a  perturbation  expansion,  such  that 
the  computationally  expensive  push  matrices  (involving  trigonometric  function  evaluations)  are 
again  computed  once  in  each  computational  cell,  while  each  particle  is  transported  according  to  a 
hybrid  scheme  involving  the  one  described  here,  and  a  rapid  scheme  such  as  the  Boris  algorithm 
for  the  small  perturbation  SB  ox  a  similar  procedure  that  does  not  involve  time  staggering.  This 
can  be  investigated  in  the  future. 

6.  References 

[1]  J.  P.  Boris,  “Relativistic  plasma  simulation  -  optimization  of  a  hybrid  code”,  Proc.  4‘^  Conf. 
Numerical.  Simulation  of  Plasmas,  pp.  3-67,  J.P.  Boris  and  R.A.  Shanny  eds.  Naval 
Research  Laboratory,  Wash.  D.C.,  Nov.  2-3  1970. 

[2]  C.  K.  Birdsall  and  A.  B.  Langdon,  Plasma  Physics  via  Computer  Simulation,  pp.  58-63,  Inst, 
of  Phys.  Publ.,  Bristol,  UK,  1991. 

[3]  W.  W.  Lee,  Phys.  Fluids  26(2),  556-562  (1983). 


14 

Distribution  A;  Public  release,  distribution  unlimited. 


Appendix  A:  Boris  Algorithm 


The  Boris  algorithm  is  defined  by  the  following  steps,  from  the  velocity  at 
fields  at  t : 


1.  v”  =v(i-At/2)  +  -^— E(0 

m  2 

2.  v'=v" +-^— v^xB(0 

m  2 


3.  v+=v“+- 


q  At 
m  2 


1  + 


q  At 
m  2 


-v'xB(0 


B 


4.  v(i+At/2)  =  v^ +-^— E(0 

m  2 

The  position  is  advanced  by  the  additional  step: 

5.  r{t+Ai)  -  r(t)  +  At  \{t+Atl2) 


t  -  At  1 2  and  the 

(AT) 

(A.2) 

(A.3) 

(A.4) 

(A.5) 


The  algorithm  could  also  be  written  a  different  way.  Let  us  define  the  following  vector: 


(A.6) 


where  co^-  —  B  is  the  cyclotron  frequency  (unsigned)  and  b  =  B  /  B  is  the  unit  vector  along 
m 

the  magnetic  field.  Steps  2  and  3  can  be  combined  into  the  form: 


2  _  2 

v  =v - ^Pxv  + - ^pxpxv 


(with  p  =  |p| ).  The  equivalent  matrix  form  is: 


(A.7) 


V  =  V  +■ 


\+p^ 

or  equivalently: 

... 


"  0 

+A 

-Py^ 

'-Pl-Pl 

PxPy 

PxPz  ' 

-A 

0 

+Px 

•V  + - - 

\+p^ 

PyPx 

-Pl-Pl 

PyPz 

[+A  “A 

0  7 

,  PzPx 

PzPy 

-Pl-Pl  ^ 

-V  +- 


\+p^ 


'  0 

+Pz 

-Pyl 

0 

'PxPx 

PxPy 

PxPz" 

-Pz 

+Py 

0 

-Px 

+A 

0  7 

_  A 

•  V  + 

\+p^ 

PyPx 

yPzPx 

PyPy 

PzPy 

PyPz 

PzPz, 

(A.8) 


(A.9) 


The  origin  of  the  Boris  algorithm  is  made  clear  by  the  following.  Consider  the  rotation  step  as 
follows: 


leading  to: 


or  in  matrix  form: 


Jv  =  v^-v-  (v-+Av)xB(0 

m  2 


+  q  At  At 

V  +A - Bxv  =v  -A - Bxv 


m  2 


m  2 


(A.10) 

(A.11) 
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r  1 

-Pz 

+Py'" 

r  1 

'^Pz 

-Py\ 

^Pz 

1 

-Px 

•V+  = 

-Pz 

1 

+Px 

VPy 

^Px 

1 

J 

^+Py 

-Px 

1 

(A.12) 


N_ 

The  matrix  on  the  LHS  can  be  inverted  to  yield: 

^  1  +  ^x  '^Pz^PxPy  ~Py^PxP: 
~Pz'^PvPx  1  +  Pi:  '^Px'^PvP: 


n:'=-  ^ 


\+p^ 


y+Py+PzPx  -Px^PzPy 


1  +  A 


yr'z 

2 


(A.13) 


The  product  •  N_  is: 


n;'-n_=-  ^ 


\+p^ 

which  can  be  decomposed  into  the  form: 


\+Pl-Pl-Pl  ^^Pz^^PxPy  -^Py^^PxPz 
-iP^^iPyP^  \+Pl-Pl-Pl 
+2^+2^,^,  -IP^^IPJ,  \+Pl-Pl-P 


(A.  14) 


'y  J 


n:'  •n_=^(i)+  — 


\+p^ 


\+p^ 


0 


^Pz  ~Py 
-Pz  0 
y^Py  -Px  0 


+  ■ 


l+y0^ 


PxPx 

PxPy 

PxPz 

PyPx 

PyPy 

PyPz 

PzPx 

PzPy 

PzPz 

(A.  15) 


We  see  that  this  is  equivalent  to  (A.9),  and  therefore  the  steps  2  and  3  of  the  Boris  algorithm  are 
equivalent  to  a  time-centered  scheme  for  the  gyro-motion  (A.  10). 


Note  that  the  Boris  algorithm  is  operator-splitting  the  electric  acceleration  from  the  magnetic 
rotation.  We  could  also  look  into  a  complete  operator  definition  without  this  splitting,  by 
considering  the  full  time-centered  algorithm: 


h+1/2  ’  H-V2 


a  q  tyt 

=  — AAE„  + - -(v„+i/2+v„_^2)xB, 

m  m2 


(A.  16) 


which  becomes 

N+-v„^i/2  =a„  At  +  N_-v„_^2  (A.17) 

where  a„  =  {qlm)  ■  E„ .  The  solution  is  already  expressed  using  the  matrices  of  (A.13)  and  (A.  15). 
To  simplify  the  notation,  let  us  define  the  following: 

,  ,  r  0  -Py^ 

Mo= - y(l),  Mi= - - 


-Pz  0  +Px 


^Pv  ~Px  0 


and  M2= 


1 


l+y0" 


V, 

then  the  solution  is  expressed  as: 

v„.v2=N;‘-kAi)+(N;'N_)-v 

with: 

N;'  =Mo  -tMi  -tMj 

and  N;'N_  =(l-/?^)Mo -t2Mi -t2M2 


n-l/2 


PxPx 

PxPy 

PxPz] 

PyPx 

PyPy 

PvPz 

(A.  18) 

PzPx 

PzPy 

PzPz] 

(A.19) 

(A.20a) 

(A.20b) 

In  the  case  where  the  magnetic  field  is  aligned  along  the  z  -axis,  considerable  simplification 
occurs.  We  can  easily  see  that  (A.9)  leads  to  the  following  relations  for  the  parallel  and  transverse 
components  respectively: 

(A.21a) 


V|l  =Vm 
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and  v|  - 

where  5  -co^bd .  This  leads  to: 


1 

1+(J/2)2 


Av  =  V„^1/2-V„-1/2 


=  (-v^+ia„A0-(-v  -|a„A0 

i-d^H  5  o' 

l+(A/2)  0  0  0 

V 

ox,  keeping  terms  of  order  only: 

a^At  +  ((y^A0(v,+|a^A0  -\{co^tdf  v, 
Av  =  <  a^At  -  ((y^A0(v,+Ta^A0 
a^At 


(A.22b) 


(A.23a) 


(A.23b) 
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